A Universal Energy Functional for Trapped Fermi Gases with Large Scattering Length 
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Yoram Alhassid conjectured that the total energy of a harmonically trapped two-component Fermi 
gas with large scattering length is a linear functional of the occupation probabilities of single-particle 
energy eigenstates. We confirm his conjecture and derive the functional explicitly. We show that 
the functional applies to all smooth potentials having a minimum, not just harmonic traps. We also 
calculate the occupation probabilities of high energy states. 
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Introduction— It is well known that the energy of a 
noninteracting Fermi gas in a trap can be expressed as a 
summation over single-particle energy eigenstates 



E = T.^. 



(1) 



where v and a label the orbital and spin states respec- 
tively, en is the single-particle energy level, and n,j„ is 
the occupation probability of the state {u, a). 

For interacting systems, however, there is no general 
relation between the total energy, which includes inter- 
action energy, and n^a alone. 

In this Letter we study a remarkable exception to this 
general rule, in a Fermi gas with strong interactions. 

We consider a two-component (a ='f, \.) Fermi gas in 
which the s-wave scattering length a and other relevant 
length scales, such as the mean inter-particle spacing d 
and the thermal de Broglie wave length, are all much 
larger than the range of the interaction r^ — >■ 0. Ap- 
plications of this s-wave contact interaction model range 
from ultracold thin vapors of neutral atoms near broad 
Feshbach resonances to the neutron gas with density and 
temperature below nuclear scales, or from extremely cold 
(10"^-10"'^K) to extremely hot (lO^K) matter. 

The energy of such a Fermi gas was found to be [iHsl 



E 



h^I 



(fik 



Airma ' J (2tt)^ 






2m 



(Pk- - I) 



(2) 



where pko- is the momentum distribution, ^^(r) is the 
spatial density distribution, m is each fermion's mass, 
V{r) is the external potential, and I = limfc_j.oo k'^Pka 
is the contact, a parameter characterizing the number 
of small pairs of fermions [l|, |2[. The contact is at the 
center of many universal relations for fermions with s- 
wave contact interaction |ll434|. 

Equation ^, containing two continuous distributions 
[pkcr and ncr{r)], is considerably more complicated than 
Eq. ([T]) which involves merely the occupation of single- 
particle levels. In the case of a harmonic trap, ni,a- are 
just a discrete set of numbers. 



Yoram Alhassid conjectured that for harmonically 
trapped fermions with large scattering length, the total 
energy might still be a linear functional of n^a j35| . 

We will show that his conjecture is true by deriving 
the functional explicitly. We also show that it is valid 
for all smooth potentials having a minimum (which we 
set to zero without loss of generality) . One can equally 
well apply this functional to anharmonic traps, periodic 
potentials (eg, an optical lattice for cold atoms), or dis- 
ordered potentials, etc. 

This universal energy functional is 



E = 
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where n^ = '^a'^uiy- The contact I is contained in the 
asymptotic behavior of 



Pcr(e) = ^n^cr(5(e 



at high energy, 

P<T(e) 



hi 



-3/2 



coarse grained 



0{e 



-5/2N 



(4) 



(5) 



We also calculate the occupation probabilities of in- 
dividual single-particle energy eigenstates at e^, 3> 
max{|£'|/A^, h^/ma^, h^/md^, AF}, where N is the 
number of fermions, and AV^ is the characteristic range 
of potentials involved in the A^-body state. (For exam- 
ple, if many fermions form a cloud in a trap, AV^ is the 
change of V{v) from the trap minimum to the edge of the 
cloud.) The result is 

+ 0{eZ% (6) 

where k^, = y/2rruZ/h, C(r) is the contact density [l|, |2| 
[related to the total contact: J C{r)d'^r — I], (j),^{r) is 
the normalized wave function of the I'-th single particle 
orbital, satisfying the Schrodinger equation 



[ - ihy2m)V^ + T/(r)] M^) = e.'^.(r), 



(7) 



D(r) is the "contact current" [see Eq. dTT]) below], and 
j,(r) = {h/m)lm^lV^, ~ 0[(e,/m)i/2|<^,|2] is the 
probability current of the ly-th single-particle orbital 
state. 

Equations ([3]), ([5]) and ([S]) apply to both energy eigen- 
states and thermal ensembles, both equilibrium and non- 
equilibrium states, both few-body and many-body sys- 
tems, both strong (|a| > d) and weak (|a| <C d) in- 
teractions, both symmetric (iV-|- = iV^) and polarized 
{Nf ^ Ni) states. In the many-body regime, they 
are valid for all phases, including normal and superfluid 
phases. 

In the following we first derive an expansion for the 
one-particle density matrix Po-(r,r -I- b) = {%pl{r)il;„{r + 
b)) at a small separation b. From this we derive Eq. ([3]) 
(exploiting the propagator of a single particle in a short 
imaginary time), Eq. (O, and Eq. ([6]). The derivations 
are for energy eigenstates but can be easily extended to 
thermal ensembles and non-equilibrium states. 

One-Particle Density Matrix Consider a normal- 
ized iV-body energy eigenstate {N — N^ + N^) 

!</.) =(7Vt! m^'^^ J DlDi(f>{n, • • • , r^^. Si, • • • , s^J 
X i,l{r,) ■ ■ ■ ^|(rAr^)^|(si) • • • VI(sArJ |0) , (8) 

where |0) is the particle vacuum, V'J(i") is the standard 
fermion creation operator, and we have introduced short- 
hand notations dJ = U^d^r^ and DJ = U%d^s^. 
When ri and Si are close, 0(r i , • • • , r^v^ , Si , • • • , s^^ ) sat- 
isfies the Bethe-Peierls boundary condition 

(p^A{^{ri -|-Si);r2---rjv^S2---SArJ(|ri -Si|"^ ^ «"^) 
+ 0(|ri-si|). (9) 

We now expand the 1-particle density matrix 

Pt(r,r + b) = (0|^j(r)V't(r + b)|0) 

^N^ dId\ (j)* (r, r2 • • • rw^si • • • sjvJ0(r + b, ra • • • ) 

(10) 

through order 0{b^) at a small distance b. Because of the 
singularity of when two fermions in different spin states 
are close (see above), we divide the 3{N^ + Ni — 1) dimen- 
sional integration domain into region TZ^ (in which every 
spin-down fermion lies outside of the sphere of radius e 
centered at r, namely |s^ — r| > e ior fi — 1,- ■ ■ , N^) and 
its complement, TZ^. Here e is small but e> b. In TZe we 
expand 0(r-f b,r2 • • -rAr^si • • • sn^) in powers of b, while 
in 7?.e we use Eq. ([9]) which is sufficient for evaluating the 
integral in 7?,e through order b^ . In TZ^ it is possible for 
two or more spin-down fermions to come inside the small 
sphere of radius e centered at r, but the contributions 
from such cases are suppressed by Fermi statistics and 
are of higher order than 0{b^). When the integrals in 



the two regions TZ^ and TZ^ are added, all dependencies 
on e are canceled, yielding the following clean expansion: 

Pt(r, r -F b) = nt(r) + C(r)(-6/87r -I- b^ l2ATTa) + b • u^(r) 
- 37r6b • w(r)/2 - nbh ■ w*(r)/2 

3 

+ Y.v^,,{v)hb,/2 + 0{b^), (11) 

where bi is the i-th Cartesian component of b, 

«tW = ^T I DlD\\^{rv2 ■ ■ ■ TN.Si ■ ■ ■ s^Jp (12) 

is the density of spin-up fermions at position r, 

C(r) = 16n^ N^N^ J DlDi\A{r; r2 ■ ■ ■rN,S2- ■ ■ sjvJP 

(13) 

is the contact density [ij, [^ , 

w(r) = N^Ni f dId^^ A* {r;r2 • • • r^^S2 • • • satJ 

X VrA(r;r2 •••rAr^S2---SAr^) (14) 
is related to the center-of-mass motion of small pairs, and 

ut(r) EE 7V^ lim / DlDi(j)*irr2---rN,Si---SN^) 

xS/r(j)irr2---rN^Si---SN^), (15) 

wtij(r) = N^ lim / D|D];0*(rr2 • • • tat Si • • • sn,) 

(rr2---rAr^Si---SAr^). (16) 



dndrj 



There is of course a completely analogous expansion for 
P4,(r,r -I- b) involving the same C(r) and w(r). 

In addition to the contact density C(r) [l|, I4I, we in- 
troduce a ^^ contact currenf^ 



D(r) = (87r^n/m)Imvir(r) 



(17) 



There is in general no continuity relation between the 
contact density and the contact current, because the 
small pairs may dissociate or associate. 

Universal Energy Fiinctional— For any {N-f + Ni)- 
body energy eigenstate |0) and any (3 satisfying Re/3 > 
we define an absolutely convergent series: 

Since the fermion annihilation operator 

C. = /ciV<^:(r)V.(r), (19) 



we have 



Ja{P)= / d'rd'r'Up{Y,v')p„{r,r'), 



(20) 



where C/^(r, r') = ^^ e~^^''0y(r)(/)* (r') is the propagator 
of a single particle moving in the potential ^(r) within a 
time —ih/3. For a small positive /?, at |r — r'| ^ h^/pjrn 
the propagator is exponentially suppressed, while at 
|r — r'l ^ hy/pjm we have a "short imaginary-time ex- 
pansion" 

C/^(r,r') = (27rft2/3/™)-3/2{i _ [y(i.) + y(r')]/3/2} 

X exp [ - m(r - r')V2?i^/3] + 0(^^/^), (21) 

provided that V{t) is smooth. But when |r — r'| is small 
we also have a systematic expansion for P(j(r,r') [see 
above]. Substituting both expansions into Eq. pO|) we 
obtain a systematic expansion for Ja{P) at small /3: 



J^(/3) =N„~ (fiZ/47r2)v/27r/3/TO + h^ip/^Trma 
-13 f (frV{r)n„{r) 



3 

+ ife(3/2m)j dV^«,,,(r)+0(/33/2) 



(22) 



From the A^-body Schrodinger equation 

/J— 1 m'— 1 

= £;(/), if r^ 7^ Sp/ for aU ^, ^', 
one can show that 



2m 



(23) 



^ /dV[nr)n.(r)-^E^-^«l =^' (24) 



so the summation of Eq. 



over (7 yields 






Applying ^ to the above expansion, defining p{t) 
J2va "i^i^a-Sie — €1,), and taking /? — ?► 0, we find 



(25) 
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(26) 



and thus Eq. ([3]). 

Asymptotics of Pa{()— When /3 == it/h is purely 
imaginary and the real "time" t — > ±0, the net contribu- 
tion to the integral in Eq. (|20p from |r — r'| ^ ^JW(t\Jrn 
is exponentially small because of the rapid oscillation of 



the propagator U. When |r — r'| ~ yjh\t\/m the ex- 
pansion in Eq. ([2T|) with j3 replaced by it/h holds [36[. 
Therefore, Ja{it/h) has a "short real-time expansion" by 
simply setting /3 — it/h in Eq. ([22l) . and thus has an 
0(1/111) singularity at t = 0. So the function P(j(e) de- 
fined in Eq. Q, which is the Fourier transform of J^: 



Pa{€) = {2nhy 



Ja{it/h)e"''/^dt, (27) 



has a "coarse-grained" asymptotic formula at high en- 
ergy shown in Eq. ([5]). For a deep trap, the fact that 
Pa(() remains a discrete sum of delta functions at large 
e, rather than turning into a continuous curve, can be 
traced to the singularities of J„{it/h) at nonzero i's. But 
in a "coarse grained" distribution, /9o-(e)| coarse grained = 
/_ 9(e')Pcr{e ~ e')de' , where the convolution factor may 
be chosen as g{e') — exp(— e'2/A^)/(A-\/7r) with a large 
width A, these singularities are "washed out", because 
Per (e) I coarse grained IS thc Fouricr transform of the prod- 
uct oig{t) and Ja{it/h), where g{t) is the inverse Fourier 
transform of g{e) and decays exponentially at |t| ^ h/X. 
For the validity of Eq. ([5]) , the energy resolution A should 
not grow faster than constant x y/e. 

Asymptotics of n^o-^ From Eq. (fT9|) we find 



dV0,(r) dH(f>:ir + h)p,{r,r + h). (28) 



At large e,y, the integrand as a function of b oscillates 
rapidly. The only significant contribution comes from 
the power-law singularities of Po- at b — ^ 0. According 
to Eq. (fTT|) . the leading order singular term is ex |b|, for 
which we write 0*(r + h) = {~h'^/2me^fVt(j)l{r + b) 
with relative error ^ 0(e~^) according to Eq. ([7]) [37|. 
Integration by parts over b yields ex / d^b<j)^{r + h)5{h), 
leading to the first term on the right hand side of Eq. ([6]) . 

The next order singular term in p^ is ex 6b, for which 
we write (/)*(r-hb) = {-h^ /2me^fVl(j)l{v + h) with rel- 
ative error ^ 0(e~^). Integration by parts over b yields 
ex j d^bcj)l[v + h)Vb5{h) = -V0*(r). Further integrat- 
ing by parts over r, omitting contributions ~ 0{e~^), and 
using Eq. ^ and the identities Rew(r) = VC(r)/327r2 
and V • jiy(r) = 0, we obtain the second term on the right 
hand side of Eq. ([6]) . 

Because any single-particle orbital state and its time 
reversal have the same energy but opposite probability 
currents, the second term on the right hand side of Eq. ^ 
has no net contribution to the distribution Pai^)- 

We now illustrate Eq. ([6]) with a symmetric unitary 
(jfeiT-aj ^ 1) Fermi gas at zero temperature, confined by a 
spherical harmonic trap of angular frequency to. At large 
N the local density approximation (LDA) for the contact 
density is valid: C{r) = Cikp^^l — r^/R'^)'^. Here kpo is 
the local Fermi wave number at the trap center, R is the 
LDA cloud radius, and Cikp is the contact density of the 
homogeneous unitary Fermi gas with Fermi wave number 



kp. For a high energy orbital ly = {j,l,mz) with energy 
^i' = (j + ^/2)hu!, orbital angular momentum quantum 
number I, and magnetic quantum number m^ we find [38| 



16Ci7V5/6 



Stt 31/6^3/4^5/2 
U/3 



1- 



/(? + l) 



^V^my/^j 



5/2 



(29) 



if Z(; + l) < Ay^{3N)^'^j. For higher Z, the classical orbit 
is outside of the LDA cloud radius, and n^a becomes ex- 
ponentially small. Here ^ is the Bertsch parameter [39J, ie 
the ratio between the ground state energy of the unitary 
Fermi gas and that of the noninteracting Fermi gas at 
the same density. In latest numerical and experimental 
studies, ^ < 0.38 and Ci « 0.12 [si-lallill- 

To conclude, we have shown that the total energy of 
fermions with large scattering length (|a| ^ re) in any 
smooth potential is a simple linear functional of the occu- 
pation probabilities of single-particle energy eigenstates. 
We have also derived asymptotic expressions for the occu- 
pation probabilities of high energy states. These results 
can be verified experimentally by measuring the energy 
and the occupation probabilities independently. They 
also provide robust constraints on theories of trapped 
Fermi gases, including fermionic atoms in optical lattices. 
These results can be extended to lower dimensions, to 
fermions with unequal masses, and to bosons and Bose- 
Fermi mixtures. 

The universal energy functional Eq. ([3]) implies a new 
approach to the difficult many-body problem at large 
scattering length: by identifying nontrivial constraints 
on the occupation probabilities, one can minimize the 
functional to find the many-body ground state energies 
in external potentials. 
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